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We present a parameter estimation procedure based on a Bayesian framework by applying a 
Markov Chain Monte Carlo algorithm to the calibration of the dynamical parameters of a space 
based gravitational wave detector. The method is based on the Metropolis-Hastings algorithm and 
a two-stage annealing treatment in order to ensure an effective exploration of the parameter space 
at the beginning of the chain. We compare two versions of the algorithm with an application to a 
LISA Pathfinder data analysis problem. The two algorithms share the same heating strategy but 
^Nj . with one moving in coordinate directions using proposals from a multivariate Gaussian distribution, 

while the other uses the natural logarithm of some parameters and proposes jumps in the eigen- 
space of the Fisher Information matrix. The algorithm proposing jumps in the eigen-space of the 
Fisher Information matrix demonstrates a higher acceptance rate and a slightly better convergence 
towards the equilibrium parameter distributions in the application to LISA Pathfinder data . For 
q " this experiment, we return parameter values that are all within ~ la of the injected values. When 

y-^ | we analyse the accuracy of our parameter estimation in terms of the effect they have on the force- 

per-unit test mass noise estimate, we find that the induced errors are three orders of magnitude less 
f*"^ ■ than the expected experimental uncertainty in the power spectral density. 
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I. INTRODUCTION 



LISA Pathfinder [H-U is a European Space Agency mission aiming to develop and test the key technology for the 
forthcoming space-based gravitational wave detectors [3| • A considerable amount of the LISA Pathfinder mission time 
will be dedicated to a series of experiments for the calibration of the model parameters for the spacecraft dynamics. 
Such procedures will be periodically replicated on future space-based gravitational wave detectors in order to ensure 
that the sensitivity of the instrument stays at the values required for the measurement of gravitational waves. Since 
the instrument calibration has a direct impact on the scientific outcome of the mission, it is of fundamental importance 
£ — ' to develop a powerful and solid procedure for parameter estimation that allows an accurate calibration in a short 
amount of time. 

The problem of the calibration of parameters describing the dynamics was already discussed in previous works 
\ focusing on LISA Pathfinder 0-0] . In Ref H, Q a frequentist approach to the problem was developed using both 
t-H ■ linear and non-linear parameter estimation techniques respectively. In Ref Q a Bayesian framework was developed for 
the second LISA Pathfinder mock data challenge, based on a classic Metropolis-Hastings Markov Chain Monte Carlo 
. i-H . (MCMC) scheme for parameter estimation 0]. From the point of view of the present paper, the method presented 
in Ref [3] is comparable with the algorithm working with the direct physical parameters and proposing the MCMC 
?H . jumps from a multivariate Gaussian distribution. 

The MCMC algorithm is a stochastic method that is very useful for carrying out Bayesian inference. The algorithm 
allows us to easily sample from a sometimes unknown probability distribution. By generating an ergodic Markov chain, 
we can sample from the underlying target distribution. The longer the MCMC, the more statistically independent 
the chain is and the better we model this distribution. In general MCMC algorithms are relatively easy to implement. 
While guaranteed to converge to the target distribution, the convergence time of the MCMC is problem dependent 
and is virtually impossible to estimate. However, the convergence can be accelerated if we use proposal distributions 
for the chain which we believe to be close to the target distribution. MCMC algorithms are normally divided into 
two parts : a "burn-in" phase where we allow the chain to find an equilibrium state, and a true MCMC which we use 
afterwards for statistical purposes. The end of the burn-in phase is usually decided "by-eye" as it is problem specific. 

We propose a Bayesian approach to the calibration of the parameters of the spacecraft dynamics that is based 
on Metropolis-Hastings MCMC. The Metropolis-Hastings MCMC is a reference method for parameter estimation in 
Bayesian Statistics (see for example Ref [10|]) and in space-based Gravitational waves data analysis (for an overview, we 
refer the reader to Ref [ll[ and references therein) where it has demonstrated its flexibility and power in the estimation 
of parameters for a variety of GW sources [l2[ ■ We compare here two different algorithms for the construction of the 
Markov Chain. One works with the physical parameters and proposes MCMC jumps from a multivariate Gaussian 
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distribution, the other uses a conversion from some of the physical parameters to their natural logarithm in order 
to calculate the Fisher Information Matrix. This improves the numerical stability for the inversion operation to 
obtain the covariance matrix, and proposes MCMC jumps in the eigen-space of the Fisher Information matrix. Both 
methods implement a double-stage annealing heat treatment in order to effectively explore the parameter space at 
the beginning of the chain. 

The rest of the paper is laid out as follows. In section|TI]we provide details on the LISA Pathfinder dynamics and on 
the experiment for the calibration of parameters of the dynamics along the principal measurement axis. In section ITTT1 
we develop the mathematical basis of our algorithm, while in section IIVI we report the results of the application 
of our methods to LISA Pathfinder and we analyze the two methods in terms of acceptance rate, convergence to 
the equilibrium distribution and error on the test mass noise estimation caused by the uncertainty on the estimated 
parameters. 



II. A TEST CASE FOR LISA PATHFINDER MISSION 



LISA Pathfinder is a controlled three-body system composed of two test masses and the enclosing spacecraft. One 
test mass is free-falling along the principal measurement axis and it is used as a reference for the drag-free controller 
of the spacecraft. The second test mass is actuated at very low frequencies (below I mHz) in order to follow the 
free falling test mass. This actuation scheme provides a measurement bandwidth of 0.5 < / < 100 mHz in which 
both the test masses can be considered to be effectively free-falling. The system has two output channels along the 
principal measurement axis, which sense the displacement of the spacecraft relative to the free-falling test mass and 
the relative displacement between the individual test masses. We focus our attention on the principal measurement 
axis dynamics for which we can build a model that contains at least the five key parameters reported in Table U 6] . 



TABLE I: Parameter descriptions and the simulation and parameter range values used in the simulation and in the MCMC 
search. 'Simulation values' are the values used to produce the set of synthetic data, they represent the 'true' values for the 
parameters in our test. The 'search range' defines the extent of the uniform prior distribution used in the MCMC. 



Name 



Description 



S&i 



Actuation gain of the drag-free control 
loop acting on spacecraft thrustcrs 

Actuation gain of the electrostatic sus- 
pension actuating the second test mass 
at low frequency 

Interferometer cross-talk coefficient cou- 
pling the first channel with the differen- 
tial channel 

Stiffness coupling the motion of the 
spacecraft to the motion of the first test 
mass 

Stiffness coupling the motion of the 
spacecraft to the motion of the second 
test mass 



Simulation Value 



Search Range 



1.05 
1.05 



[0.5, 1.5] 
[0.5, 1.5] 



1 X 10" 



[-1 x 10~ 3 , 1 x 10 3 



1.4 x 10~ D s" 



[1 x 10 , 1 x 10~ 5 ] s" 



2.2 x 10~° s" 



[1 x 10~ 7 , 1 X 10~ 5 ] s 



We simulated two classical experiments [6| for LISA Pathfinder: one stimulating the motion of the spacecraft with 
respect to the free-falling test mass and the other forcing a movement of the actuated test mass with respect to the 
free-falling mass. The signal in the first experiment is obtained by forcing the drag-free controller of the spacecraft 
to change sinusoidally its set-point. This sends a command to the thrusters that produces a sinusoidal movement 
of the spacecraft with respect to the free-falling test mass. In the second experiment we force a sinusoidal bias in 
the set-point of the system that controls the low frequency actuation of the second test mass. As a consequence a 
sinusoidal voltage is applied by the actuators to the second test mass and a movement with respect to the free-falling 
test mass is induced. For the generation of our synthetic signals we adopted the standard scheme in which the output 
of the instrument can be written as s(t) = h(t) + n(t), where h(t) is the time-dependent noiseless output of the 
LISA Pathfinder model for the dynamics (the model is described in Ref @) and n(t) is an additive noise contribution 
generated in accordance to the procedure described in Ref [f3j]. The noise n(t) is a Gaussian distributed correlated 
noise that reproduces the current knowledge of the LISA Pathfinder noise budget [l4j]. This noise is principally 
composed of thruster noise in the channel sensing the displacement of the spacecraft relative to the free falling test 
mass (first channel) and of capacitive actuation noise on the channel sensing the relative displacement between the 
test masses (differential channel). These two channels produce a displacement signal at the output that needs to be 
converted into force-per-unit-mass for the analysis of the noise affecting the test masses. The conversion is based on 
a model of the spacecraft dynamics that is written in terms of the parameters described in Table HI The calibration of 
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FIG. 1: 2D slices in the 5D likelihood surface for the current parameter set. Note that some of the axes are in log-scale while 
others are in a linear scale. We clearly see that over the parameter ranges, the surface contains a single global solution. 



the parameters of such a model therefore has a direct impact on the LISA Pathfinder test-mass noise budget, which 
is the main scientific target of the LISA Pathfinder mission. 

In order to better understand our parameter estimation problem, we calculated the 5D likelihood surface over a 
uniform grid of points spanning the parameters ranges defined in Table [I] 2D slices of the likelihood surface are 
presented in Figure ([1} where we clearly see that there is one global solution over the 'search range'. The 'search 
range' for the parameters is defined from the prior knowledge of the system that principally comes from construction 
constraints and requirements. 



III. BAYESIAN ANALYSIS AND MARKOV CHAIN MONTE CARLO 



Bayesian analysis has become a very useful tool in the field of gravitational waves [], and has also been applied to 
LISA Pathfinder data analysis Q • One of the advantages of Bayesian analysis is that it provides a direct method of 
calculating the marginalised posterior density function (PDF) for each of the parameters via Bayes' theorem 

- n{\)p(s\X) 
p{s) 

where p(X\s) is the posterior density function (PDF) for the solution A given the data s(t), 7r(A) denotes the prior 
probability of the parameters and p(s\X) is the likelihood, which we will define below. The constant in the numerator 
is the Bayesian evidence and is defined by 

P (S) = [ d\7T(\) P ( S \\), (2) 



which as can be seen in independent of the solution set A and is unimportant for our study. 
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While previous works have used sophisticated variants of the MCMC to search over a wide parameter space, in 
this study, as we can see from Figure ([1]), we are dealing with a simplified surface. Even starting at the boundaries 
of our parameter priors, we can see that the surface contains a unimodal solution. Because of this we can use a 
straightforward MCMC. Starting with a detector output s(t) = h(t) + n(t), where hit) is the induced signal and n(t) 
is the noise in the detector, we generate a template h(t; Xi) by choosing a random starting point in the parameter 
space Xi, and then draw from a proposal distribution in order to propose a jump to another point in the space Ai+i. 
To compare both points, we evaluate the Metropolis-Hastings ratio 

H = 7r(Ai + i)p(s|Ai + i)g(Ai|Ai + i) 
TT{X i )p{s\X l )q(X l+ i\X l ) 

where we now define p{s\X) as the likelihood by 

p(s\X) = C (X) = Ce-H*!-*, (4) 

where C is a normalization constant. We arc unconcerned by this constant as we ultimately arc interested in the 
likelihood ratio between the original and proposed points, hence cancelling the constant. The quantity q{X\y) is the 
transition kernel or proposal distribution used for jumping from A^ to A; + i. Once a jump is proposed, it is then 
accepted with probability a = mm(l, H), otherwise the chain stays at A^. While we are sure that the MCMC will 
eventually converge to the global solution, the timescale of this convergence is unclear. So, while we can never be 
sure that our starting values are close to the true values, it can take a long time for the chain to "walk" to the global 
maximum. A common way of decreasing this travel time is to heat the likelihood surface using a technique called 
simulated annealing. This has the effect of lowering and fattening structures on the surface. We will describe this 
method in more detail at a later stage. 

For the LISA Pathfinder experiments, we expect to have an extremely high signal-to-noise ratio (SNR) (~ 10 5 — 10 6 ) . 
In this case the detector output has the form s(i) w h(t). For this high SNR limit, it is expected that the errors in 
the parameter estimation will have a Gaussian probability distribution given by 

where V ^ is the Fisher information matrix (FIM) 



r — 



- / dh 



and r = det{T^ v ). The angular brackets above denote the noise- weighted inner product for two real functions s(t) 
and hit) 

( h = 2 r -A* \kf)~s*(f) + h*(mf)] , (7) 



o S n (f) 
where 

h(f)= / dth(t)e 2mft (8) 



is the Fourier transform of the time domain waveform hit) and an asterisk denotes a complex conjugate. The quantity 
S n (f) is the one-sided noise spectral density of the detector. For the LISA Pathfinder project, there is no closed-form 
solution to the Fisher matrix calculation as it is not possible to take the derivatives of the template with respect 
to the parameters, nor evaluate the integral analytically. Therefore, the derivatives in the inner product need to be 
calculate numerically using a one-sided difference algorithm. Also, as the templates are expensive to generate (on the 
order of seconds), we update the FIM calculation every 100 iterations during the annealing scheme and as we believe 
that the FIM should remain almost constant at the global solution, every 5000 iterations during the MCMC. 

While we have seen that the likelihood surface is highly simplified compared to other gravitational wave problems, 
in that there exists a single mode solution, our goal here is to present an algorithm that can be started from any point 
in the parameter space. In order to accelerate the convergence of the chain to the global peak, we use a mixture of 
annealing methods as presented in [l5|, [l6[ . This first means re- writing the likelihood as 



(X) = C exp (-± (s - h (X) | s _ h (a) , (9) 



5 



where we introduce a heat scaling j3. The first annealing scheme used is called thermostated annealing. While the 
SNR for the experiments are large, and even though there is no chance of getting stuck on a secondary solution in 
this particular problem, it is difficult to choose the initial annealing temperature by hand. If the temperature is too 
high, all structures on the likelihood surface are completely flattened and the MCMC aimlessly explores a flat surface. 
Therefore, as it is again costly to generate templates, we do not want to waste computer cycles by choosing the initial 
temperature too high. On the other hand, in the case of a multimodal solution, if the temperature is too low, the 
chain can get stuck on a secondary maximum and never reach the global solution. Therefore, we allow the algorithm 
to chose its own initial annealing temperature according to 




< SNR < SNRt 



SNR 1 



(10) 



SNR T 



SNR > SNRt 



which ensures that once we attain a SNR of greater than a threshold SNR, SNRt, the effective SNR never exceeds 
this value. In our algorithm we run the thermostated annealing for the first 2000 iterations. For this problem, we 
choose our threshold SNR as follows. As previously stated, the SNRs involved in this problem are very large. If we 
take the boundary values for the parameters, as given in TablcU the SNR at the boundary is 1.473 x 10 6 , while at the 
injected values we have a SNR of 1.486 x 10 6 . To choose the threshold SNR, we arbitrarily demand an initial value of 
/3 = 500. Solving for the threshold SNR in the above expression then provides a threshold value of SNRt = 6.58 x 10 5 . 
After this we use a second annealing scheme, which is standard simulated annealing, to cool the surface down. To do 
this we use an annealing of the form 

f lO^ 1- *) < i < T c 

(11) 

where i represents the chain index, and T c represents the number of chain points during the cooling schedule. For this 
study, £ = 2o<7io(/3rocL) i s ^ ne initial simulated annealing heat- index defined by the maximum temperature at the end 
of the thermostated annealing phase 0^_ x - Again, as there is no risk of getting stuck on secondary maxima, we cool 
down the surface much quicker that we normally would for other problems. Here we set the cooling schedule equal to 
T c = 3000 iterations. 

One of the things that we wanted to test in this study was the difference between moving in eigen-directions and 
coordinate directions when we update the parameters. In Ref [HI, EH the proposed updates were done in the following 
manner. Once the Fisher matrix is calculated, we find the eigenvalues and eigenvectors V^„. We then update the 
parameters in the eigen-directions according to 

/ l~a\ 5 t/ 

(12) 





where D is the dimensionality of the problem and n v — A/"„(0, 1) is a vector of normal distributed random numbers 
with zero mean and unit variance. The factor of 1/vD ensures that the total jump in the parameter space has a 
length of ~ 1 a, while the heat factor f3 ensures that while the temperature is high, we take larger steps in the chain 
and thus faster explore the parameter space 

To move in coordinate directions, we first invert the Fisher matrix to obtain the variance-covariance matrix, i.e. 

c = c^ = (r^r 1 . (13) 

As C is a real symmetric matrix, we then use a Cholesky decomposition to write the covariance matrix as C = LL*, 
where L* is the Hermitian conjugate or adjoint matrix of L, and L is a unique lower triangular matrix. We should 
note here that for a real symmetric matrix the Hermitian conjugate is equivalent to the matrix transpose. We then 
generate a vector of zero mean, unit variance random Gaussian values n = n„, and we update the parameter values 
according to 



A" -> A^ + y/0 (L.n) , (14) 

where the update is along the coordinate directions for the parameters. 

In the analysis conducted in this article, we combine the information provided by two LISA Pathfinder output 
channels and two separate experiments. The inner product in Equation |7|) can be generalized in a straightforward 
way to the case of multiple output channels by 
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FIG. 2: The instantaneous acceptance rate for two Markov chains moving in eigen-directions (blue-solid) and coordinate 
directions (orange-dashed). We observe a factor of two improvement by moving in eigen-directions, as opposed to moving in 
coordinate directions. 



s-h(A) s-h(A)) =2 / S 1 -h 1 (X),...,S Nc -h Nc (X) S~ (/) 
o 



si - fei(A) 

SNc — /liVc(A). 



df + cc (15) 



Here £„(/) is the cross-spectral matrix, which contains the power spectral density at frequency / on the diagonal 
and the cross-spectral density on the off-diagonal terms, Nc is the number of output channels. Since we assume that 
the outputs of two distinct experiments arc independent, the total likelihood for a number iV e of experiments is just 
the sum of the likelihoods from each experiment, i.e. 



Ci 



N, 



(16) 



IV. RESULTS 



Our study involves two different types of MCMC. For the movement in eigen-directions, we choose to work with the 
parameter set A M = {Adf,Aif s ,S&\,\n wfjln By working in the natural logarithm of the parameters {u;f,a;§} 

we lower the condition number of the Fisher matrix, making it numerically more stable for inversion purposes. This 
matrix is still not as stable as we would like due to the small values of the parameter Sai- However, as this parameter 
can also have negative values we cannot use its logarithm. For the coordinate direction chains, we use the physical 
parameter values. Each chain was run for 2 x 10 5 iterations, with 5,000 iterations of annealing. To confirm the relative 
simplicity of the problem, starting with injected values of {1.885, 2.0399, -0.00746, 9.135 x 10~ 6 , 6.327 x 10~ 6 }, both 
chains converged to the global maximum within 3000 iterations. However, just to be safe, we will label the first 8,000 
chain points as "burn-in" . These points will be rejected when it comes to constructing a statistical analysis from the 
chains. These initial values correspond to a separation of approximately {10 3 , 10 5 , 10 4 , 10 4 , 10 4 }<r from the injected 
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FIG. 3: The instantaneous parameter means for two Markov chains moving in eigen-directions (blue) and coordinate directions 
(orange). The two methods show similar rate of convergence and despite we are dealing with an unimodal problem it takes 10° 
samples to have a reasonable convergence. 



values, where a here is the standard deviation calculated from the chain and presented later in Tabic HT1 In Figure ([2]) 
we plot the instantaneous acceptance rates for both the eigen-direction chain (blue-solid) and the coordinate direction 
chain (orange-dashed). We can see that we gain roughly a factor of two by using jumps in eigen-directions rather 
than coordinate directions with the chains having final acceptance rates of 63% and 31% respectively. This means 
that to independently sample the same posterior density, the coordinate direction chain would need to be twice as 
long as the eigen-direction chain. 

In Figures (j3|) and (j4|) we plot the convergence of the means and standard deviations of each parameter from the 
two MCMC methods. We can see from comparing both methods that they have a similar rate of convergence. We 
should point out here that we would not expect both chains to find exactly the same value. If we were to run iV 
chains, we would expect the chains to produce ./V values all roughly within ±ler of the true value. Our criterion of 
evaluation here is how quickly the chains flatten with no large fluctuations. Although after 2 x 10 5 iterations, it is 
still difficult to prefer one method over the other. This is because of the fact that even though we are dealing with 
a simplified problem, it still takes 10 5 iterations for the Metropolis-Hastings algorithm to begin to converge. This 
suggests that for full convergence, the chains would need to be run for longer. However, it is still clear that we are 
finding the global peak very quickly. 

In Figure (0 we plot the marginalised, normalised histograms from the MCMC. On the x-axis we plot the true 
value offset scale (A# — AeO/ocr, where \r is the true injected parameter value, Ac is the chain value and oc is the 
standard deviation as calculated from the chain. We can see that both chains give us distributions very close to the 
true values. If we look closer, we can however see that the coordinate direction distributions have a higher kurtosis 
due to the lower acceptance rate. In all cases, the distributions are almost Gaussian with a skewness of ~ 10~ 2 . 
While there may seem to be no reason to prefer one method over the other, we will focus on the results from the 
eigen-direction chain due to its higher acceptance rate. We present the results in full in Table HU but we quote here 
that the chain returned mean values for each parameter that were within {1.17, —0.58, 0.44, —0.23, — 0.69}cr of the 
true answer. 

One of the aspects of the problem that we can also investigate from the MCMC is the correlations between 
parameters. Using the chains, we calculate the variance-covariance matrix using 



C^v — 



1 



A'* 



(17) 



where N = 192, 000 is the number of chain points (ignoring the 8,000 burn-in points at the beginning of the 2 x 10 D 
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FIG. 4: The instantaneous parameter standard deviations for two Markov chains moving in eigen-directions (blue) and coor- 
dinate directions (orange). 




FIG. 5: This Figure displays the marginalised ID histograms from the chains moving in eigen-directions (left) and coordinate 
directions (right). The x-axis origin is set at the true value of the parameters and the offset is the scaled difference between 
the true value Xr and the chain value Ac, divided by the standard deviation of the chain ac- Both chains give distributions 
very close to the true values but from a closer observation the coordinate direction distributions have a higher kurtosis due to 
the lower acceptance rate. 



iterations) . An overbar denotes the mean of the parameter as calculated from the chain. Once we have the variance- 
covariance matrix, we can calculate the matrix of standard deviations and correlation coefficients according to 



/! = V 



(18) 



where <j b 



'C w are the standard deviations of the marginalised chains for each parameter. The numerical values 
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TABLE II: Simulated and recovered parameter values from the eigen-direction MCMC. The final column gives the standard 
deviation scaled offset from the true value. 

Simulation Value Mean Standard Deviation (\r — \c)/&c 



1.05 


1.049445 


4.7356 x 10~ 


-4 


1.17 


1.05 


1.050002 


3.3727 x 10- 


-6 


-0.58 


10" 5 


9.9314 x 10~ 6 


1.5669 x 10 


-7 


0.44 


1.4 x 10" 6 


1.400041 x 10" 6 


1.8048 x 10" 


10 


-0.23 


2.2 x 10" 6 


2.200111 x 10" 6 


1.5974 x 10" 


10 


-0.69 



TABLE III: Test mass noise response to parameters variations. The response of the noise as a function of each parameter is 
linear n{v) = Av + B, therefore we obtained the slope A by a linear fit to the data. The noise variation An at la is obtained 
multiplying A by the sample standard deviation of the parameters obtained from the chain. 



A 


A 


An at la [ms" 2 Hz" 1/2 ] 


Ad f 


1.9 x 10" 28 


9 x 10" 32 


Aifs 


6.9 x 10" 15 


2 x 10" 20 


Sai 


3.5 x 10~ 16 


5 x 10~ 23 




2 x 10~ 12 


4 x 10~ 22 




1.2 x 10" 11 


2 x 10~ 21 



for the matrix, where the order of the rows and columns correspond to {A^f, Aif s , Sai, ^1} 



are 



D„ 



( 4.736 x 1CT 4 -9.921 x 1CT 3 1.662 x 1CT 3 -6.891 x 1CT 3 
3.373 x 1CT 6 1.259 x 1CT 2 6.697 x 10" 1 
1.567 x 1CT 7 -2.267 x 1CT 1 
1.805 x 1(T 10 

V 



-1.547 x lO- 3 ^ 
7.659 x 1CT 1 
8.097 x 10~ 3 
8.789 x 10- 1 
1.597 x 10~ 10 / 



(19) 



We can see that the parameter set is mostly uncorrelated. However, we do see strong correlations between Aif s and 
the two parameters (w 2 ,w 2 ), i.e. 0.67 and 0.77, and between w 2 and themselves, i.e. 0.88. 

It is interesting here to quantify the accuracy of our result in terms of the LISA Pathfinder test mass noise budget. 
As discussed above, the parameters of the dynamics play an important role in the analysis of the test mass noise, since 
they are used for the calibration of the model that converts the displacement noise to force-per-unit-mass. On the 
basis of such a model we calculated the noise response to the parameters as reported in Table Hill This data is obtained 
calculating the noise-per-unit-mass on an uniform grid of values for the parameters. The grid is centered around the 
simulation 'true value' and extends over 10cr, where a is the sample standard deviation obtained by the chain. The 
response of the noise as a function of each parameter is linear n(v) — Av + B, therefore we obtained the slope A 
reported in Table Hill by a linear fit to the data. The noise change An at la is obtained by multiplying A by the values 
of <7 obtained by the chain. As can be seen in Table ITTll even combining the results for all the five parameters we obtain 



A? 



2xl0~ 21 [ ms" 2 Hz 



-l/2i 



Even assuming that we use 6 months data (the full mission duration) this value is three 



orders of magnitude lower than the expected experimental uncertainty (2.3 x 10~ 18 [ms^Hz^ 1 / 2 ]). For this calculation 
we calculated the power spectral density with the well known Welch's overlapped segment averaging method, using 
a 4 sample Blackman-Harris window [17| and 50% segment overlapping. A reasonably long LISA Pathfinder noise 
run can last for 3 days, and the current estimation for the noise level is 7.2 x 10~ 15 [ms~ 2 Hz -1 / 2 ] at 1 mHz [Tij ]. 



_1/2 ]. While 
' 1 ' 2 ], where 

the uncertainty is simply calculated as the value for the power spectral density divided by the number of averaging 
segments. 



Using this information in calculating the expected uncertainty, we obtain a value of 1.4 x 10 [ ms Hz 
if we assume a 6 month long data run, we obtain an experimental uncertainty of 2.3 x 10 -18 [ms _2 Hz 



V. CONCLUSIONS 



We have presented a study of two different MCMC methods for an experiment with LISA Pathfinder. While both 
methods are based on the Metropolis-Hastings algorithm, and both use the Fisher matrix to generate the covariance 
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matrix, they differ in their treatment of the parameters and the way in which they move in the parameter space. 
In one case, we work with the physical parameter values and make MCMC jumps in the coordinate directions. For 
the other chain, we generate a more numerically stable Fisher matrix by working in the logarithm of two of the five 
parameters. This lowers the condition number of the matrix and makes it more stable during the inversion process to 
obtain the covariance matrix. We demonstrate that while the experiment chosen for this study is quite simple, in that 
a singe global solution exists in the likelihood surface, we gain roughly a factor of two in acceptance rate by moving 
in eigen-directions rather than in coordinate directions. Our chains both end up with parameter estimations that are 
very close to the true values. This is to be expected due to the extremely high SNR of the signal. However, as we 
are conducting a Bayesian inference, we are mostly interested in the full posterior distribution for each parameter. In 
this case, due to the higher acceptance rate, it is clear that calculating a more stable Fisher matrix and moving in 
eigen-directions is an optimal choice. The accuracy of the recovered parameters has been analyzed in terms of the 
effect they have on the force-per-unit-mass test mass noise estimate. We demonstrated that the error introduced by 
the parameters estimated from the MCMC chain is in any case more than three orders of magnitude less than the 
expected experimental uncertainty in the power spectral density. 
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